/////////////////////////////////////////////////////////////////////////////
*Figure 3: Child earnings penalty
/////////////////////////////////////////////////////////////////////////////
//Note: This code is adapted directly from the replication code provided by Kleven, Landais, and Sogaard (2018)



clear all


gen gender = "`gender'"
save "temp_eventstudy.dta", replace

//Indicate window for event study
local start -3
local end 5
local base -2


set more off
foreach gender in women men{
	foreach var in earnings extensive intensive wagerate regular1 {

	
		display "    " 
		display "`gender' `var' `group'"
		display "   " 
		
//Child penalty
//Keep only those with at least one first child born during the sample, and a window of at least 2 years before and after
use Data\klips_workingdata_hh, clear
drop if years<2000
bysort pid: egen firstchildborn = total(hh_firstchildborn)
drop if firstchildborn==0

//Keep selected gender
keep if gender=="`gender'"

//Generate event time tracker
gen eventyear = years if hh_firstchildborn==1
bysort pid: egen year_firstchildborn = mean(eventyear)
gen eventtime = years-year_firstchildborn
local totalyears = `end'-`start'+1
drop if eventtime<`start'|eventtime>`end'

//Generate balanced panel
bysort pid: egen mineventtime = min(eventtime)
bysort pid: egen maxeventtime = max(eventtime)
drop if mineventtime>`start' | maxeventtime <`end'

//Create dummies
char eventtime[omit] `base'
xi i.eventtime


		//Generate dependent variable for regression
		if "`var'"=="earnings"{
			gen var = earnings
		}
		if "`var'"=="extensive"{
			gen var = working
		}
		if "`var'"=="intensive"{
			gen var = worktime 
		}
		if "`var'"=="wagerate"{
			gen var = wage 
		}
		if "`var'"=="regular1"{
			gen var = regular1
		}

		
		//Run regression 
			reg var _Ieventtime* i.years i.age [aw=_lweight98_], robust cluster(pid)
	

		predict var_p, xb

		gen b  = .
		gen bL = .
		gen bH = .
		replace b  = 0                                         if eventtime == `base'
		replace bL = 0                                         if eventtime == `base'
		replace bH = 0                                         if eventtime == `base'

		local c = -`start'+1
		local newbase = 1 + `base' - `start'
		
		foreach i of numlist 1(1)`totalyears' {
			display "Eventtime " `i'-`c'
			if `i' != `newbase' {
				replace b  = _b[_Ieventtime_`i']                             if eventtime == `i'-`c'
				replace bL = _b[_Ieventtime_`i'] - 1.96*_se[_Ieventtime_`i'] if eventtime == `i'-`c'
				replace bH = _b[_Ieventtime_`i'] + 1.96*_se[_Ieventtime_`i'] if eventtime == `i'-`c'
				}
			}
		

********************************************************************************
* CREATING COUNTERFACTUAL
********************************************************************************
		gen var_c  = var_p - b
		gen var_cL = var_p - bL
		gen var_cH = var_p - bH
		
		keep eventtime var_* b bL bH
		sort eventtime
		collapse var* b bL bH, by(eventtime)

		gen variable = "`var'"
		gen gender   = "`gender'"
		gen group = "`group'"
		
		append using "temp_eventstudy.dta"

		save "temp_eventstudy.dta", replace

}
	}



********************************************************************************
* FIGURE 3 
********************************************************************************



foreach var in earnings extensive intensive wagerate regular1 {
	display "    " 
	display "`var'"
	display "   " 

	use "temp_eventstudy.dta", clear

	keep if variable == "`var'"

		gen gap    = (var_p - var_c) /var_c
		gen boundL = (var_p - var_cL)/var_c
		gen boundH = (var_p - var_cH)/var_c

	keep b var_c gap boundL boundH variable eventtime gender

	reshape wide b var_c gap boundL boundH, i(variable eventtime) j(gender, string)
	
	local commonopts xlabel(`start'(1)`end', nogrid) `yaxis' xline(-.5)  ysc(r(-0.7 0.55)) ylabel(-0.60(0.20)0.40)
	
	if "`var'"=="earnings"{
		global label Earnings
		global ylabel Earnings
		global axis1 `commonopts'
	}
	if "`var'"=="extensive"{
		global label Participation Rate
		global ylabel Participation Rate
		global axis1 `commonopts'
	}
	if "`var'"=="intensive"{
		global label Hours Worked
		global ylabel Hours Worked
		global axis1 `commonopts'
	}
	if "`var'"=="wagerate"{
		global label Hourly Wage
		global ylabel Hourly Wage
		global axis1 `commonopts'
	}
	if "`var'"=="regular1"{
		global label Regular Job
		global ylabel Regular Job
		global axis1 `commonopts'
	}	

	twoway (rarea boundLmen   boundHmen   eventtime, color(gs14%50) lcolor(gs12%50) lwidth(vvthin))        ///
		   (rarea boundLwomen boundHwomen eventtime, color(gs14%50) lcolor(gs12%50) lwidth(vvthin))        ///
		   (connected gapmen   eventtime, lcolor(gs8) mcolor(gs8))      			    			 ///
		   (connected gapwomen eventtime, lcolor(black) mcolor(black))  ///
	, graphregion(color(white)) xtitle("Event Time (Years)", size(small))  ///
	legend(order(3 "Male" 4 "Female") col(2) ring(0) position(6) size(small)) ///
	$axis1 title("$ylabel", size(medsmall))

	graph save Childevent_`var'.gph, replace

}

grc1leg Childevent_earnings.gph Childevent_intensive.gph Childevent_extensive.gph  Childevent_wagerate.gph, cols(2) legend(Childevent_earnings.gph) note("All graphs show variable percent difference in the effect of children relative to the counterfactual of no children." "We estimate the impact relative to two years before arrival of first child)")
graph export Figures\Figure3.pdf, replace

graph combine Childevent_regular1.gph
graph export Figures\AppendixFigure1.pdf, replace

/*
//Erase intermediate files
erase temp_eventstudy.dta
foreach file of local erasefiles{
	erase `file'
}
